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Abstract. In this work, wc present a numerical method to consistently ap- 
proximate solutions of a spatially discrete, double sine-Gordon chain which 
considers the presence of external damping. In addition to the finite-difference 
scheme employed to approximate the solution of the difference-differential 
equations of the model under investigation, our method provides positivity- 
preserving schemes to approximate the local and the total energy of the system, 
in such way that the discrete rate of change of the total energy with respect 
to time provides a consistent approximation of the corresponding continuous 
rate of change. Simulations are performed, first of all, to assess the validity 
of the computational technique against known qualitative solutions of coupled 
sine-Gordon and coupled double sine-Gordon chains. Secondly, the method 
is used in the investigation of the phenomenon of nonlinear transmission of 
energy in double sine-Gordon systems; the qualitative effects of the damping 
coefficient on the occurrence of the nonlinear process of supratransmission are 
briefly determined in this work, too. 



1. Introduction 

The well-known sine-Gordon equation is a partial differential equation tliat ap- 
pears in many applications, either in its original form or as a slight modification 
of the classical version. For instance, a damped sine-Gordon equation appears in 
the study of long Josephson junctions between superconductors when dissipative 
effects are taken into account [1]. A similar partial differential equation with dif- 
ferent nonlinear term appears in the study of fluxons in Josephson transmission 
lines [2]. Meanwhile, a modified Klein-Gordon equation appears in the statistical 
mechanics of nonlinear coherent structures — such as solitary waves — , in the form 
of a Langevin equation (see [3], pp. 298-309). 

The spatially discrete version of the sine-Gordon equation also has many impor- 
tant applications. For instance, a coupled system of discrete sine-Gordon equations 
may describe a chain of harmonic oscillators coupled through Hookean springs [?] 
or a system of Josephson junctions attached through superconducting wires [5]. In 
the former case, a system initially at rest, with void initial velocities and sinusoidal 
Dirichlet boundary condition at one end, is employed to study the phenomenon 
of supratransmission of energy [4j [6j , which is a nonlinear process characterized 
by a sudden increase in the amplitude of wave signals generated by the perturbed 
boundary, for driving amplitudes above a critical value called the supratransmission 
threshold. This phenomenon is also present in the investigation of the transmission 
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of energy in chains of Josephson junctions, except that, in this case, a harmonic 
Neumann boundary condition needs to be imposed upon the physical problem for 
the sake of meaningfulness [5^ . 

We must remark, in this point, that several other bounded, nonlinear regimes 
present the phenomenon of supratransmission when they are harmonically per- 
turbed at one end, such is the case of Klein-Gordon arrays [4], Fermi-Pasta- Ulam 
systems , Bragg media in the nonlinear Kerr regime [8] , and even in spatially con- 
tinuous, bounded media described by undamped sine-Gordon equations [S]- More- 
over, the presence of the phenomenon of nonlinear supratransmission is also found 
in discrete, double sine-Gordon chains [6 . However, the specialized mathematical 
literature unfortunately lacks studies to approximate the occurrence of the process 
in these systems. 

Nevertheless, in this work, we study the process of nonlinear supratransmission 
in dissipative, double sine- Gordon chains, employing a dissipation-preserving finite- 
difference scheme. Of course, there exist many analytical |10[ 111) and numerical 
[l2l fT3l [T4] techniques to approximate solutions of Klein- Gordon-like equations. 
The computational method employed in this work distinguishes from many other 
techniques available in the literature in that it consistently approximates not only 
the solution of the physical model under study, but also the local energy density, the 
total energy, and the rate of change of the energy with respect to time. Moreover, 
the nonnegative character of both the energy density and the total energy — a char- 
acteristic which is not preserved by other numerical techniques [15] — , is preserved 
by our method. As we shall see later on, these qualities are, by definition, highly 
desirable characteristics of any computational technique employed in the study of 
the phenomenon of supratransmission. 

In Section [21 we present the system of ordinary differential equations that moti- 
vates our study, together with the local energy functions associated, and the total 
energy of the system. A proposition which summarizes the expression of the de- 
rivative of the total energy with respect to time is provided in this stage, as well 
as a brief description of the process of nonlinear supratransmission, particularly, in 
double sine-Gordon systems. Section [3] introduces the numerical method employed 
to approximate the solutions of the model under investigation, the local energy 
distribution, and the total energy of the system. A subsection on the numerical 
properties of the method summarizes the properties of consistency established for 
convenience in the appendices, and another subsection presents some remarks on 
the computational implementation of our technique. Section 2] presents simulations 
of sine-Gordon and double sine-Gordon systems, obtained by means of our method. 
The former regime is employed only for validation purposes, while the simulations 
on the latter (which follow the same methodology proposed in 4 ) are aimed at 
establishing the existence of the process of nonlinear supratransmission in this sys- 
tem. Finally, this work closes with a section of concluding remarks and further 
directions of research. 



2. Preliminaries 

2.1. Physical model. Throughout this work, we let Z^v = {1, 2, . . . , — 1}, for 

every positive integer N; obviously, we will assume that A^ > 1 for the sake of 
non-triviality. Moreover, we let Zjv = "^n U {0, A^}, that is, Z^r = {0, 1,2,..., A^}. 
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Let iV be a positive integer and, for every n — 0,1,..., A^, let u„ be a real 
function of time t > 0. Moreover, let c be a positive real number, and let 7 be a 
nonnegative number. Throughout, we consider a mechanical chain of nonlinear os- 
cillators obeying the system of ordinary differential equations with initial-boundary 
conditions 

^ - SPu^ + 7^ + V'{u,,) =0, VneZ^, 
u„(0) ~ 0, V?i G Zat, 

I ^(0)-0, VneZ^, 
at 

uo{t) = 4>{t), yt > 0, 

UNit) - UN-lit) ^0, Vt > 0. 

In other words, we consider a spatially discrete, bounded system initially at rest, 
with zero initial velocities, perturbed at the left end by a function (j) which we 
will assume to be continuous, and with discrete Neumann boundary condition on 
the right end. The constant 7 is immediately identified as the external damping 
coefficient, while c is clearly the coupling coefficient between nodes. Here, the 
spatial, second-difference operator 

(2) S^^^Un = C^(Un+l - 2Un + Un-l), Vn G Z^y, 

has been employed for convenience. 

For the sake of concreteness, we will consider a driving function of the form 

(3) (l3(t) = Asm{nt), 

where the driving amplitude A and the driving frequency fl are positive real num- 
bers. Moreover, we consider a potential function of the form 

(4) ^(m) = ^ - i [2cosu + cos(2u)] , 

whence the double sine-Gordon law V'{u) = ^ [sinw + sin(2u)] readily results. In- 
deed, let c = If Aa; is relatively small (or, equivalently, c is relatively large), 
then the system of ordinary differential equations of ([1} approximates the spatially 
continuous, partial differential equation 

where L — NAx, and w is a function of space x and time t. This equation is clearly 
identified with the classical, double sine-Gordon equation with constant external 
damping. 

Of course, different potentials may give rise to other important models in mathe- 
matical physics. For instance, V{u) = 1 — cos(w) is the potential for the sine-Gordon 
regime, while V{u) = — jj-u* + corresponds to the potential of a nonlin- 
ear Klein-Gordon equation. In fact, it is important to warn the reader that the 
dissipation-preserving numerical technique presented in this work is valid not only 
for the double sine-Gordon potential, but also for any differentiable function V 
defined on the real line. 
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n — 1 n + 1 

Figure 1. Forward-difference stencii for the approximation to 
tfie partial differential equation ([1} at time tk, using the finite- 
difference scheme (ITS)) . The black circles represent known approx- 
imations to the actual solutions at times tk-i and tk, and the 
crosses denote the unknown approximations at time tk+i- 

2.2. Energy of the system. Let n £ Zat. For physical reasons, it is important 
to notice that the local energy of the nth node in the undamped system governed 
by (dl is provided by the expression 

(6) - \ 

where the spatial, first-order difference operator 5x is defined through 

(7) Sj^Un = c (u„+i - u„) , Vn e Zat. 

In these terms, the total energy E of the system ([T]) is obtained by adding 
the local energies for n S Z^-, and the potential due to the coupling in the 
boundaries of the chain. In other words, 

(8) E=Y1 ^« + \ (^-^o)' • 

The following proposition is easy to establish. 

Proposition 1. The rate of change of energy with respect to time of a system 
governed by (QP is given by 

□ 

As a corollary, the system ((T)) conserves the total energy if no external damping 
is present and, either ^ is a constant function or a void Neumann condition is 
imposed on the left end of the chain. 



dUn 

dt 



V{Un 



Vn e Zat, 
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Figure 2. Graphs of approximate solution wgo versus time, of the 
60th node of a system ([T]) of length N — 200, for c = 4, 7 = and a 
potential V{u) = 1 — cosu. The system was perturbed by means of 
(|3]) with il = 0.9, and two different amplitudes were used: A = 1.77 
(top graph) and A — 1.78 (right column). The insets depict the 
corresponding temporal evolution of the local energy of the 60th 
node. 



2.3. Nonlinear supratransmission. As observed in the literature (see [4j[6l[7]), 

the double sine- Gordon system ([IJ, as well as the nonlinear Klein- Gordon and the 
sine-Gordon chains, and the classical /3-Fermi-Pasta-Ulam systems, presents the 
phenomenon of supratransmission of energy, which is a nonlinear process charac- 
terized by a sudden increase in the amplitude of wave signals propagated into a 
nonlinear medium by a driving source which irradiates at a frequency in the for- 
bidden band-gap. The mechanism of this transmission of energy is through the 
generation of localized, nonlinear modes at the driving boundary, in the form of 
moving breathers or soliton solutions 16 . 

More concretely, consider a nonlinear system of any of the types mentioned in the 
previous paragraph, which is perturbed at one end by a harmonic function of the 
form ([3|) , with a fixed value in the forbidden band-gap of the system. Relatively 
small driving amplitudes A result in the propagation of practically no energy into 
the system; however, as the value of A is increased, the existence of a critical value 
As, above which the system begins to absorb great amounts of energy from the 
boundary, is immediately noticed. 

The value As introduced in the previous paragraph, is called the supratransmis- 
sion threshold, and its existence has been analytically proved for discrete HI and 
continuous [9. sine- Gordon chains, as well as for systems of anharmonic oscillators 
[7]. However, as it was mentioned before, the study of the double sine- Gordon 
regime has been left practically unexplored. 

For our particular study, a simple analysis of the undamped, linearized system 
of differential equations in ([T]) shows that the linear dispersion relation is given by 

(10) a;2(fc) = l + 2c2(l-cosfc). 

In our simulations, the driving frequency fl will take on values in the forbidden 
band-gap region il < 1. 
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Figure 3. Graphs of approximate solution wgo versus time, of the 
60th node of a system ([T]) of length N — 200, for c = 4, 7 = and a 
potential given by Q . The system was perturbed by means of ([3j 
with fl = 0.9, and two different amplitudes were used: A = 1.03 
(top graph) and A = 1.04 (right column). The insets depict the 
corresponding temporal evolution of the local energy of the 60th 
node. 



3. Numerical method 

3.1. Finite-difference scheme. Let iV be a positive integer, and let T be a posi- 
tive real number. In order to approximate the solutions of the system ([1]) at time T, 
we fix a regular partition of the interval [0, T] of the form = to < ti < . . . < tj^j = 
T, of norm At = T/M. Additionally, we let be the numerical approximation 
of the actual value of w„ at time tk, for k = 0, 1, ... , M. Moreover, in order to 
simplify our notation, we define the temporal differences 

(11) 



(12) 6, 



u 



At 

fe+i 



(13) <5r^< 



2At 
{Atf 



for every n G Zjv and k G Z^/ . Furthermore, for such values of n and fc, we employ 
the temporal average operator 

(14) = + 

and the discrete derivative of V with respect to u and the time average of V at u, 
respectively: 

(15) ^ 



(16) 

Finally, let 

(17) 0fc = 0(tfe) 
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Figure 4. Graphs of approximate solution versus node site n and 
time t, of a system ([IJ of length N — 200, for c = 4, 7 = and 
a potential given by dH). The system was perturbed by means of 
(jS]), with J7 = 0.9, and two different amplitude values: A — 1.03 
(left) and A = 1.04 (right). 



With these conventions at hand, the numerical method to approximate solutions of 
([T]) is summarized as follows: 

(<5f ) - + 7*^^ + Si^'^v) «) = 0, Vn e Z^, 

Vn £ Zat, 
Vn £ Zat, 
VA: G Zm, 
0, ykeZM- 

For convenience, the forward-difference stencil of this method has been depicted in 

Fig.m 



(18) 




N-l 



3.2. Energy scheme. With the same notation as in the previous paragraph, the 
local energy of the system ([1]) at the nth node and at the kth time step will be 
approximated by means of the discrete formula 

(19) = I 

where n G Z^v and k e Za/. Meanwhile, the total energy of the system at time tk 
is calculated through the expression 

(20) ^'== E^n + ^E^^- 

n^Zj^ l — k 

Before closing this stage of our investigation, it is important to point out that 
the local energy function Hn in ([6]) is nonnegative for the case of the double sine- 
Gordon regime; in addition, its discrete counterpart, namely, Eq. (I19p . is likewise 
nonnegative. It follows that the total energy of the system ([T]) as given by ([S]), as 
well as the discrete total energy (PU)) are both nonnegative at any time. 



n k+1 Is 

^ ^ 

j=n—l l=k 
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Figure 5. Graphs of approximate local energy versus node site n 
and time t, of a system ([Ij of length N — 200, for c = 4, 7 = 
and a potential given by ([4]). The system was perturbed by means 
of ([3|), with r2 = 0.9, and two different amplitude values: A = 1.03 
(top graph) and A= 1.04 (right column). 



3.3. Numerical properties. As mentioned previously, the numerical method pre- 
scribed by the expressions ^TE\\ . (|19p and (PO)) preserves the positivity character of 
the local and the total energy of the system ([1]). Moreover, the finite-difference 
schemes presented in p8)) provide consistent solutions of ([T]) of order the second or- 
der in time (see Appendix |X] for a brief discussion of the consistency of the method). 

The fact that the local energy estimate (fT9|) is a consistent approximation of 
the continuous local energy (jH), and that the discrete total energy ([20l) . in turn, is 
a consistent estimation of the corresponding continuous expression is evident. 
The following result shows that this consistency is also preserved on the grounds 
of the rate of change of energy with respect to time. 

Proposition 2. Consider the finite- difference scheme \18\) . with local energy given 
by I119\) . and total energy I180\) . Then, the discrete rate of change of energy of the 
method at time tk-i is given by 

(21) S.E'^^^ = -7 E {^t'^-nY 



Proof. See Appendix IB] 



□ 



3.4. Computational remarks. Clearly, the finite-difference scheme (IT51) is non- 
linear and implicit when V is not a constant function, as it is the case of the double 
sine-Gordon chain. Thus, in order to approximate the solution of the system ([l]) at 
time tk+i when the approximations at times tk and tk-i are at hand, we employ 
Newton's method for nonlinear systems of equations. 

Once, again, let us adopt the conventions of Section [XT] For every k e Zm, 



let u*"' = (uq, u'l 
equation of p^ . that is, let 



'■N 



), and let /„ be the left-hand side of the nth difference 



(22) 
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Figure 6. Graph of approximate total energy over the time pe- 
riod [0, 200] of the undamped system ([T]) versus driving amphtude, 
subject to harmonic driving of the form ([3]) and a potential 
The parameters c — i, N — 200 and fl = 0.9 were employed in the 
simulations. 



for every n 6 Zjy. Additionally, let 



(23) 
(24) 



Moreover, let f — (/o, /i, . . . , /at). Using a recursive process, assume that the 
vectors u'"' and u''^^ have been previously computed. Then 

(25) u'^+i = u*^ - y, 

where y is the {N + l)-dimensional vector which satisfies the matrix equation 

(26) J(u^> = -f(u''). 

Evidently, the (iV + 1) x {N + 1) matrix J is the Jacobian matrix of f , which is 
given by 



(27) 



J(u^) 



/ 1 

a di a 

a d2 a 



V 



\ 





a d^-i a 

-1 1 y 
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Figure 7. Graph of approximate total energy over the time period 
[0, 200] of the undamped system ([1]) versus driving amplitude and 
driving frequency (left), subject to harmonic driving of the form 
([3]) and a potential (jlj. The parameters c = 4 and N = 200 were 
employed in the simulations. The left graph is the checkboard plot 
of the top one. 



where 

(28) a = -| 

(29) dn = T-^+c' 



(At)2 +^ +2At 



(ut'-utT 

for every n G Zjy. The tridiagonal system (PB)) is solved then employing Grout's 
reduction technique with pivoting [17) . 

Of course, for our simulations, Newton's method requires of a stopping criterion 
in order to approximate the vector u*^"*"^ in Eq. (|25l) . Particularly, in this work, 
this criterion is given by the condition ||y|j2 < e, where the tolerance parameter e 
is equal to 1 x 10~^, and || • ||2 is the classical Euclidean norm in E^+^. 

4. Simulations 

Throughout this section, we consider a system governed by ([ij, where the driving 
function assumes the sinusoidal form ([3]). In order to avoid the creation of shock 
waves produced by the sudden movement of the driving boundary at the initial 
time, we linearly increase the driving amplitude from zero to its actual value A 
during a finite period of time Tq. Particularly, in the simulations performed in this 
study, we fix Tq — 50. 

4.1. Sine-Gordon chain. As a means to verify the validity of our method, we 
consider, first of all, a discrete chain of harmonic oscillators coupled through iden- 
tical springs, in which case, the governing equations are given by ([1]), with V(u) = 
1 — cosM. Moreover, assume that the system under study is undamped, let c = 4, 
N = 200 and fix a driving frequency of 0.9. According to [3], the supratransmission 
threshold of the system occurs around the critical value = 1.78. 
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Figure 8. Graph of approximate driving amplitude As above 
which supratransmission occurs in the undamped system ([1]) ver- 
sus driving frequency f2, subject to harmonic driving of the form 
^ and a potential The parameters c = 4 and — 200 were 
employed in the simulations. 

From a computational perspective, we let A< = 0.05, and compute approxi- 
mations to the actual solution of the initial-boundary- value probme ([T]) and the 
corresponding local energy, over a time interval of length T = 200. Under these 
circumstances. Fig. [5] presents the temporal evolution of the solution and the local 
energy of the 60th node of the system for two values of the driving amplitude: 
A = 1.77 (top graph) and A = 1.78 (below graph). The results show a drastic 
change in the qualitative behavior of solutions around the proposed critical ampli- 
tude As- These results are clearly in agreement with [4], and they are considered 
as evidence in favor of both the validity of our method and the existence of supra- 
transmission in the sine-Gordon chain. 

4.2. Double sine-Gordon chain. As mentioned before, the study of the phenom- 
enon of nonlinear supratransmission of energy in the double sine- Gordon chain is 
a topic of interest that has been left aside. In this section, however, we proceed to 
compute bifurcation diagrams similar those constructed to predict the process of 
supratransmission in discrete sine-Gordon and Klein-Gordon systems (see [11 [6]). 
So, as in the previous stage of our investigation, we consider an undamped system 
governed by H]), with parameters c = A, N — 200, T — 200, ft = 0.9, and potential 
given by (H]). Computationally, let At = 0.05. 

With these considerations. Fig. [3] presents the time-dependent graphs of the 
solution and the local energy of the 60th node of the system, for two different 
values of the driving amplitude, namely, A = 1.03 (top graph) and A = 1.04 
(bottom graph). As in the case of the discrete sine-Gordon system, we observe a 
drastic qualitative difference in the behavior of the solution and the local energy of 



12 



J. E. MACIAS-DfAZ 



X 10'* 




1 1.2 1.4 1.6 1.8 2 

Amplitude 

Figure 9. Graph of approximate total energy over the time pe- 
riod [0, 200] of the system ^ versus driving amphtude, subject to 
harmonic driving of the form ([3]) and a potential ([4]). The parame- 
ters c — 4, N — 200 and = 0.9 were employed in the simulations, 
with four different values of the damping coefScient: 7 = (solid), 
7 = 0.01 (dashed), 7 = 0.02 (dash-dotted), 7 = 0.03 (dotted). 

the 60th node, around the critical value As = 1.04. Indeed, this observation is in 
agreement with the available literature [6j. 

Fig. Ulpresents the solution of the system studied in the previous paragraph, with 
respect to node site n and time t, where t e [0,200]. The graphs of the solutions 
clearly change drastically for the two driving amplitudes considered: for A = 1.03, 
the boundary obviously does not propagate wave signals into the system; on the 
contrary, the graph corresponding to A — 1.04 shows transmission of energy into the 
medium. This observation is verified in Fig. [SJ which presents the corresponding 
graphs of local energy for the two amplitudes considered. Evidently, the qualitative 
observations done in the domain of the solutions carry over to the domain of the 
local energy of the system. 

For the next step in our discussion, we define the total energy of the system ([1]) 
over the time interval [0,T] as 

(30) Et^ [ E{t)dt, 

Jq 

where E{t) is the total energy of the system at time t, given by expression ([8]). 
Clearly, Et is consistently approximated by means of the formula 

(31) = E 
where each E'' is given by (pO|) . 
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With this notation, Fig. [6] presents the total energy over the temporal interval 
[0, 200] , of a system ^ with the same parameters as above, when the driving ampli- 
tude takes on values in the interval [0.8, 1.3]. The graph evidently shows the drastic 
change in the behavior of the total energy of the system before and after the ampli- 
tude value 1.04. With this strong evidence of the existence of supratransmission in 
the double sine-Gordon chain, the critical value Ag — 1.04 is immediately identified 
as the nonlinear supratransmission threshold for = 0.9. Obviously, these results 
are in perfect agreement with [5]. 

We have performed similar simulations for several values of n in the interval 
[0.2, 1], and values of A in [0, 4], and we have obtained qualitatively identical results. 
Indeed, Fig. [7] summarizes our findings, in the form of a graph of total energy over 
the period of time [0, 200], versus driving frequency and driving amplitude. Clearly, 
for every such driving frequency, there exists a smallest driving amplitude As above 
which the system begins to absorb energy from the boundary. From here, a graph of 
As versus driving frequency is obtained and presented as Fig. [8l It is worth noticing 
that the results are in qualitative agreement with those obtained for sine-Gordon 
chains [HE]. 

Finally, it must be mentioned that the numerical method employed in this work 
is also useful in order to establish the effects of damping in a discrete double sine- 
Gordon chain governed by ((T|). Indeed, consider a system consisting of = 200 
nodes coupled through ([T]), with c = 4, harmonically perturbed by the driving 
function ((3]) with 51 = 0.8, over an interval of time [0,200], where the potential 
function V is given by (jl]). Fig. |9] presents the effect of the driving amplitude A 
on the total energy of the system, for values of A in the interval [1,2], and three 
different values of the damping coefficient, namely, 7 = 0, 7 = 0.01, 7 = 0.02 and 
7 = 0.03. The results show the expected decrease in the total energy of the system 
as 7 is increased and, moreover, they show that the supratransmission threshold is 
slightly delayed with the presence of damping. 

5. Conclusions 

In this work, we have employed a numerical method in the study of the occur- 
rence of the process of nonlinear supratransmission in a discrete chain of oscilla- 
tors coupled with identical springs. The method proposed is consistent of order 
0{{At)'^), and it is associated to a discrete scheme to approximate the local energy 
of the chain, as well as a scheme for the total energy of the system. Both energy 
schemes consistently approximate their continuous counterparts, and the method 
has the property that the discrete rate of change of energy also approximates the 
corresponding continuous rate of change. 

The method was qualitatively tested against known approximations to the oc- 
currence of the phenomenon of nonlinear supratransmission in discrete sine-Gordon 
and double sine- Gordon chains. The simulations obtained with our method are in- 
deed in excellent agreement with the results available in the literature. Moreover, 
the method was employed in the construction of a bifurcation diagram of smallest 
driving amplitude at which supratransmission starts in the undamped system, ver- 
sus driving frequency. The graph is actually in qualitative agreement with those 
found in the literature for discrete sine-Gordon and Klein-Gordon chains, which 
are systems with the same forbidden band-gap region. Moreover, when damping is 
present, our simulations show that the process of supratransmission is still present 
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in the system under investigation, and that the appearance of the critical amphtude 
value is delayed as the damping coefficient increases. 

Of course, many avenues of research still remain open. Thus, from a practical 
point of view, it is important to provide applications of the results presented in this 
work. More concretely, following [9] , it is interesting to propose applications of the 
process of nonlinear supratransmission to the design of amplifiers of weak signals, 
or to the fabrication of detectors of ultra weak pulses, as it has been done for the 
Klein- Gordon equation [51 
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Appendix A. Consistency study 

A brief consistency analysis of the finite-difference schemes (|18p reveals that the 
numerical method proposed in this work is consistent of order 0{{At)'^). In fact, 
observe that 

(33) . 

for every n £ Z^r and every k £ Zm- Moreover, 

(34) f^i''si^^u'^^6^M.+ ^-^e'^- 

Appendix B. Energy consistency 

For the sake of simplification, we introduce the following notation, for every 
n eIin and k G 1>m: 



(35) = 77 ( 



(36) J" 



2 v^n+l + '^ii) 1 
1 [OxUq) 



4 ^ 2 

l=k 



(37) h^r. = ^ E E 



j—n—l l—k 



Clearly, the term is identified with the independent term (the term which is not 
prescribed by the summation over all n G Zjv) to the right-hand side of Eq. (f20| . 



Proof of Proposition [H It is easy to check that the following identities are valid for 
every n e Z^r and k G Zjv/: 

(38) \ {5,u'^)' ~ \ {5,ut-'f = Ai, 
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(39) ^,tV{ui) - ^^tV{ut') = {si^^Viut)) {Si'^'^) At. 

It is a tedious algebraic task (though relatively easy) to verify that the following 
equalities hold, for every n G Zj^i and k e Z^: 

(40) +1 U'^ut+A U'^S^u':^ 



(41) 



As a consequence, 

E ^^^t' - E { [^i"' - ^^'^'s'i' + €'y] {<)■ 

- -7 E (^*^"-^)' 

Moreover, 

(42) = \ (b^^h^ut) [^^^H^ut) . 

The conclusion of Proposition [2] is now reached by computing StE'^^^ and simpli- 
fying. □ 
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